function ImageBinary=PlotSynergy_L145_limb_2(L1,L4,L5,epoch,FN_i,CameraAngleCapability,DU,VU,TU)

options=odeset('RelTol',1e-13,'AbsTol',1e-13);
R_sun = 696340;
h_apex = 10000;
beta = asind(R_sun/DU);
alpha = acosd(R_sun/(R_sun+h_apex));
Theta_min = 90-alpha-beta;
Theta_max = 90+alpha-beta;
%% Compute L4

S_i = L4.S_i;

temp_time_month = month(L4.epoch);

ind = find(temp_time_month==epoch);

% plot3(S(1,:),S(2,:),S(3,:))
[row,col,~] = size(FN_i);
visibility_i=zeros(row,col);
for ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if angle_temp < CameraAngleCapability
                visibility_i(ii,jj) = 1;
            end
            if round(angle_temp)==0
                L4Location=[ii,jj];
            end
        end
    end
end


ImageBinary_L4_on_disk = mat2gray(visibility_i);

S_i = L4.S_i;

temp_time_month = month(L4.epoch);

ind = find(temp_time_month==epoch);

[row,col,~] = size(FN_i);
visibility_i=zeros(row,col);
for ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if and((Theta_min<angle_temp),(angle_temp<Theta_max))
                visibility_i(ii,jj) = 1;
            end
        end
    end
end
ImageBinary_L4_limb = mat2gray(visibility_i);

% figure;
% hold on;
% mesh(ImageBinary_L4*100,'EdgeColor','none','FaceColor','flat');
% xlabel('Longitude')
% ylabel('Latitude')
% titlestring = sprintf('L4 Month = %d',epoch);
% title(titlestring)
% set(gca,'fontsize',16)
% xticks(linspace(0,row,5));
% xticklabels({'180W','90W','0','90E','180E'});
% yticks(linspace(0,col,5))
% yticklabels({'90S','45S','0','45N','90N'});
% xlim([0,row]);
% ylim([0,col]);

%% Compute L5

S_i = L5.S_i;

temp_time_month = month(L5.epoch);

ind = find(temp_time_month==epoch);

% plot3(S(1,:),S(2,:),S(3,:))
[row,col,~] = size(FN_i);
visibility_i=zeros(row,col);
for ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if angle_temp < CameraAngleCapability
                visibility_i(ii,jj) = 1;
            end
            if round(angle_temp)==0
                L5Location=[ii,jj];
            end
        end
    end
end


ImageBinary_L5_on_disk = mat2gray(visibility_i);

S_i = L5.S_i;

temp_time_month = month(L5.epoch);

ind = find(temp_time_month==epoch);

[row,col,~] = size(FN_i);
visibility_i=zeros(row,col);
parfor ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)))
            if and((Theta_min<angle_temp),(angle_temp<Theta_max))
                visibility_i(ii,jj) = 1;
            end
        end
    end
end
ImageBinary_L5_limb = mat2gray(visibility_i);

% figure;
% hold on;
% mesh((ImageBinary_L5+ImageBinary_L5_face)*100,'EdgeColor','none','FaceColor','flat');
% xlabel('Longitude')
% ylabel('Latitude')
% title('L5')
% set(gca,'fontsize',16)
% xticks(linspace(0,row,5));
% xticklabels({'180W','90W','0','90E','180E'});
% yticks(linspace(0,col,5))
% yticklabels({'90S','45S','0','45N','90N'});
% xlim([0,row]);
% ylim([0,col])

%% Compute L1

S_i = L1.S_i;

temp_time_month = month(L1.epoch);

ind = find(temp_time_month==epoch);

% plot3(S(1,:),S(2,:),S(3,:))
[row,col,~] = size(FN_i);
visibility_i=zeros(row,col);
for ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if angle_temp < CameraAngleCapability
                visibility_i(ii,jj) = 1;
            end
            if round(angle_temp)==0
                earthLocation=[ii,jj];
            end
        end
    end
end


ImageBinary_L1_on_disk = mat2gray(visibility_i);

S_i = L1.S_i;

temp_time_month = month(L1.epoch);

ind = find(temp_time_month==epoch);

[row,col,~] = size(FN_i);
visibility_i=zeros(row,col);
parfor ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if and((Theta_min<angle_temp),(angle_temp<Theta_max))
                visibility_i(ii,jj) = 1;
            end
        end
    end
end
ImageBinary_L1_limb = mat2gray(visibility_i);
%% L4+L5

ImageBinary_L145_ondisk =1*ImageBinary_L1_on_disk+...
    2*ImageBinary_L4_on_disk+...
    4*ImageBinary_L5_on_disk;

ImageBinary_L145_limb =1*ImageBinary_L1_on_disk+100*ImageBinary_L1_limb+...
    2*ImageBinary_L4_on_disk+100*ImageBinary_L4_limb+...
    3*ImageBinary_L5_on_disk+100*ImageBinary_L5_limb;

ImageBinary_L145_circle =ImageBinary_L1_on_disk+ImageBinary_L1_limb+...
    ImageBinary_L4_on_disk+ImageBinary_L4_limb+...
    ImageBinary_L5_on_disk+ImageBinary_L5_limb;

Window = figure('OuterPosition',[100,100,900,500]);
hold on;
mesh(ImageBinary_L145_ondisk,'EdgeColor','none','FaceColor','interp');
map = [0.3 0.3 0.3 %gray 1
    0.0 0.0 1.0 % blue 2
    1.0 0.0 0.0 % red 3
    1.0 0.0 1.0 %Magenta
    0.0 1.0 0.0 % green
    0.0 1.0 1.0 % cyan
    1.0 1.0 1.0 %Magenta
    ];
colormap(map);
mesh(6*ImageBinary_L1_limb,'EdgeColor','none','FaceColor','interp');
mesh(6*ImageBinary_L4_limb,'EdgeColor','none','FaceColor','interp');
mesh(6*ImageBinary_L5_limb,'EdgeColor','none','FaceColor','interp');
xlabel('\lambda [deg]')
ylabel('\phi [deg]')
ylabel('Heliographic Latitude [deg]','Interpreter','latex')
xlabel('Heliographic Longitude [deg]','Interpreter','latex')
titlestring = sprintf('L4 Month = %d',epoch);
% title(titlestring)
set(gca,'fontsize',12)
xticks(linspace(0,row,5));
xticklabels({'180W','90W','0','90E','180E'});
yticks(linspace(0,col,5))
yticklabels({'90S','45S','0','45N','90N'});
xlim([0,row]);
ylim([0,col])
% colormap("turbo");
% cd=colorbar();
% cd.Ticks=[40,110,185,260];
% cd.TickLabels={'0','1','2','3'};
% ylabel(cd,"# of observable sats",'FontSize',15,'Rotation',270)
% cd.Label.Position(1)=4;
% plot3(earthLocation(2),earthLocation(1),1000,'ko','MarkerSize',15,'MarkerFaceColor','k')
% plot3(L4Location(2),L4Location(1),1000,'ks','MarkerSize',15,'MarkerFaceColor','k')
% plot3(L5Location(2),L5Location(1),1000,'kd','MarkerSize',15,'MarkerFaceColor','k')



% figure;
% hold on;
% [xs,ys,zs] = ellipsoid(-0,0,0,696340,696340,696340,2000);
% for ii = 1:row
%     for jj = 1:col
%         if ImageBinary_L1_limb(ii,jj) == 1
%             ImageBinary_L145_ondisk(ii,jj) = 6;
%         end
%         if ImageBinary_L4_limb(ii,jj) == 1
%             ImageBinary_L145_ondisk(ii,jj) = 6;
%         end
%         if ImageBinary_L5_limb(ii,jj) == 1
%             ImageBinary_L145_ondisk(ii,jj) = 6;
%         end
%     end
% end
% surface(xs,ys,zs,ImageBinary_L145_ondisk,'EdgeColor','none','AlphaData',0.1,'FaceAlpha',0.9);
% colormap(map);
% axis equal
% colormap(map);


end